import ps1
import re
import sys
import itertools

Seq = ps1.readSeqsFromFile("CpG-islands.txt")


    
p = 0.98; q = 0.999; nSymb = 8;
    
T = [[0.180*p, 0.274*p, 0.426*p, 0.120*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.171*p, 0.368*p, 0.274*p, 0.188*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.161*p, 0.339*p, 0.375*p, 0.125*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [0.079*p, 0.355*p, 0.384*p, 0.182*p, (1-p)/4, (1-p)/4, (1-p)/4, (1-p)/4],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.300*q, 0.205*q, 0.285*q, 0.210*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.322*q, 0.298*q, 0.078*q, 0.302*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.248*q, 0.246*q, 0.298*q, 0.208*q],
     [(1-q)/4, (1-q)/4, (1-q)/4, (1-q)/4, 0.177*q, 0.239*q, 0.292*q, 0.292*q]]
     
def argmax(array):
    max = array[0]; index = 0;
    for i in range(1,len(array)):
        if(array[i]>max): 
            max = array[i]; index = i;
        
    return index

# akl: A+ C+ G+ T+ A- C- G- T
def E(i,c):
    
    if(c=='A'):
        if(i==0 or i==4): return 1;
        else: return 0;
        
    elif(c=='C'):
        if(i==1 or i==5): return 1;
        else: return 0;
    
    elif(c=='G'):
        if(i==2 or i==6): return 1;
        else: return 0;
        
    elif(c=='T'):
        if(i==3 or i==7): return 1;
        else: return 0;

    
#E = [[1, 0, 0, 0, 1, 0, 0, 0],
#     [0, 1, 0, 0, 0, 1, 0, 0],
#     [0, 0, 1, 0, 0, 0, 1, 0],
#     [0, 0, 0, 1, 0, 0, 0, 1]]


for iseq in range(0,1):
# FORWARD PASS
    print(); print(); print(); 
    print("Gene " + str(iseq) + " : Length " + str(len(Seq[iseq])))
    print(Seq[iseq]);

    V = list(); TZ = list(); 
    V = [[0,0,0,0,0,0,0,0]];
    for k in range(0,nSymb):
        V[0][k] = E(k,Seq[iseq][0])*1/8;
    
    
    
    # FORWARD
    for i in range(0,len(Seq[iseq])-1):
        V.append([0,0,0,0,0,0,0,0]); 
        temp = [0,0,0,0,0,0,0,0];
        for l in range(0,nSymb):
            for k in range(0,nSymb):
                temp[k] = V[i][k]*T[k][l];
            V[i+1][l] = E(l,Seq[iseq][i+1])*max(temp);
        sumV = sum(V[i+1]);
        for l in range(0,nSymb):
            V[i+1][l] = V[i+1][l]/sumV ;   
        
    Z = list(); 
    for i in range(len(Seq[iseq])-1,-1,-1): Z.append(-1);
    
    Z[len(Seq[iseq])-1] = argmax(V[len(Seq[iseq])-1]);
    # BACKWARD
    for i in range(len(Seq[iseq])-2,-1,-1):
        for k in range(0,nSymb):
            temp[k] = E(Z[i+1],Seq[iseq][i+1])*V[i][k]*T[k][Z[i+1]];
        Z[i]=argmax(temp);
        
    
    print(); print();
    print("Most Probable Path: " + str(Z));
    del V; del Z;
